\documentclass[11pt]{article}
\usepackage[framed,numbered,autolinebreaks,useliterate]{mcode}
\usepackage{graphicx,booktabs,subfigure,listings,geometry,lmodern,color,xcolor,textcomp}
\usepackage{amsfonts,amsmath,siunitx,bm}
\usepackage{hyperref}
\usepackage{tikz}
\usetikzlibrary{external}
\tikzexternalize[
	figure list=true, 
	prefix=figures/
]







\geometry{left=2.5cm,right=2.5cm,top=2.5cm,bottom=2.5cm}

\title{Modeling a Bus Suspension System in Matlab and Simulink}
\author{Deng Zhongzhu}

\begin{document}

\maketitle
\tableofcontents{}


\section{Task Introduction}

\subsection{Physical Setup}

Designing an automatic suspension system for a bus turns out to be an interesting control problem. When the suspension system is designed, a 1/4 bus model (one of the four wheels), as one in Figure\ \ref{fig:suspensionSystem}, is used to simplify the problem to a one dimensional spring-damper system.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{springdampersystem.pdf}
	\caption{Model of Bus Suspension System}\label{fig:suspensionSystem}
\end{figure}

Where:
\begin{itemize}
	\item body mass \(m_1 = \SI{2500}{kg} \)
	\item suspension mass \(m_2 = \SI{320}{kg}\)
	\item spring constant of suspension system \(k_1 = \SI{80000}{N/m} \)
	\item spring constant of wheel and tire \(k_2 = \SI{500000}{N/m} \)
	\item damping constant of suspension system \(b_1 = \SI{350}{N.s/m} \)
	\item damping constant of wheel and tire \(b_2 = \SI{15020}{N.s/m} \)
	\item control force \(u\) is the force from the controller we are going to design
\end{itemize}





\subsection{Design Requirement}

A good bus suspension system should have satisfactory road holding ability, while still providing comfort when riding over bumps and holes in the road. When the bus is experiencing any road disturbance (i.e.\  pot holes, cracks, and uneven pavement), the bus body should not have large oscillations, and the oscillations should dissipate quickly. Since the distance \(x_1-w\) is very difficult to measure, and the deformation of the tire \(x_2-w\) is negligible, we will use the distance \(x_1-x_2\) instead of \(x_1-w\) as the output in our problem.

The road disturbance \(w\) in this problem will be simulated by a step input. This step could represent the bus coming out of a pothole. We want to design a feedback controller so that the output \(x_1-x_2\) has an overshoot less than \( 5\% \) and a settling time shorter than 5 seconds. For example, when the bus runs onto a \(\SI{10}{cm}\) high step, the bus body will oscillate within a range of \( \SI{+- 5}{mm}\) and return to a smooth ride within 5 seconds.



\section{Task Analyze}

By this report we are going to use three traditional analyze methods of dynamic mechanical system: \emph{deferential equation}, \emph{transfer function equation} and \emph{state space matrix} to describe the system. 

\subsection{Deferential Equation}

From the Figure\ \ref{fig:suspensionSystem} and Newton's law, we can obtain the dynamic equations:
\begin{align}
	m_1 \ddot x_1 & =-b_1(\dot x_1 -\dot x_2 )-k_1(x_1-x_2)+u                                   \\
	m_2 \ddot x_2 & =b_1(\dot x_1 -\dot x_2 )+k_1(x_1-x_2)+b_2(\dot w - \dot x_2 )+k_2(w-x_2)-u
	\label{eq:1}
\end{align}


\subsection{Transfer Function Equation}
Assume that all of the initial condition are zeros, so these equations represent the situation when the bus's wheel go up a bump. The dynamic equations above can be expressed in a form of transfer functions by taking \emph{Laplace Transform} of the above equations. Here we let \( \mathbf{X}(s)= \begin{bmatrix} X_1(s) \\ X_2(s) \\ \end{bmatrix} \), the derivation from above equations of the \emph{Transfer Functions} \(G_1(s)\) and \(G_2(s)\) of output, \(x_1-x_2\) and two inputs \(u\) and \(w\), are as follows:
\begin{align*}
	(m_1s^2+b_1s+k_1)X_1(s)-(b_1s+k_1)X_2(s)              & =U(s)                \\
	-(b_1s+k_1)X_1(s)+(m_2s^2+(b_1+b_2)s+(k_1+k_2))X_2(s) & =(b_2s+K_2)W(s)-U(s)
\end{align*}
\[
	\begin{bmatrix}
		m_1s^2+b_1s+k_1 & -(b_1s+k_1)               \\
		-(b_1s+k_1)     & m_2s^2+(b_1+b_2)s+k_1+k_2 \\
	\end{bmatrix}
	\mathbf{X}(s)=
	\begin{bmatrix}
		U(s)                \\
		(b_2s+k_2)W(s)-U(s) \\
	\end{bmatrix}
\]
\[
	\bm{A}=
	\begin{bmatrix}
		m_1s^2+b_1s+k_1 & -(b_1s+k_1)               \\
		-(b_1s+k_1)     & m_2s^2+(b_1+b_2)s+k_1+k_2 \\
	\end{bmatrix}
\]

Find the inverse of matrix \(\bm{A}\) and then multiple with inputs \(U(s)\) and \(W(s)\) on the right hand side:
\[
	\mathbf{X}(s)=\frac{1}{\det{\bm{A}}}
	\begin{bmatrix}
		m_2s^2+(b_1+b_2)s+k_1+k_2 & b_1s+k_1             \\
		b_1s+k_1                  & m_1s ^2+b_1s+k_1 \\
	\end{bmatrix}
	\begin{bmatrix}
		U(s)                \\
		(b_2s+k_2)W(s)-U(s) \\
	\end{bmatrix}
\]
\[
	\mathbf{X}(s)=\frac{1}{\det{\bm{A}}}
	\begin{bmatrix}
		m_2s^2+b_2s+k_2 & b_1b_2s^2+(b_1k_2+b_2k_1)s+k_1k_2                    \\
		-m_1s^2         & m_1b_2s^3+(m_1k_2+b_1b_2)s^2+(b_1k_2+b_2k_1)s+k_1k_2 \\
	\end{bmatrix}
	\begin{bmatrix}
		U(s) \\
		W(s) \\
	\end{bmatrix}
\]

When we want to consider input \(U(s)\) only, we set \(W(s) = 0\). Thus we get the transfer function \(G_1(s)\): \[G_1(s) =\dfrac{X_1(s)-X_2(s)}{U(s)} =\frac{(m_1+m_2)s^2+b_2s+k_2}{\det{\bm{A}}}\]

When we want to consider input \(W(s)\) only, we set \(U(s) = 0\). Thus we get the transfer function \(G_2(s)\): \[G_2(s) =\dfrac{X_1(s)-X_2(s)}{W(s)}=\frac{-m_1b_2s^3-m_1k_2s^2}{\det{\bm{A}}}\]


\subsection{State Space Matrix}
According to the differential equation in\ \ref{eq:1}, we assume 4 state variables are respectively,

\[
	\mathbf{x}=
	\left \{
	\begin{aligned}
		x_1 & = x_1                        \\
		x_2 & = x_2                        \\
		x_3 & = \dot x_1                   \\
		x_4 & = \dot x_2 -\frac{b_2}{m_2}w \\
	\end{aligned}
	\right.
\]

\[
	\mathbf{\dot x }=
	\left \{
	\begin{aligned}
		\dot x_1 & = x_3                                                                                                                                             \\
		\dot x_2 & = x_4+\frac{b_2}{m_2}w                                                                                                                            \\
		\dot x_3 & = -\frac{k_1}{m_1}x_1+\frac{k_1}{m_1}x_2-\frac{b_1}{m_1}x_3+\frac{b_1}{m_1}x_4+\frac{1}{m_1}u+\frac{b_1b_2}{m_1m_2}w                              \\
		\dot x_4 & =\frac{k_1}{m_2}x_1-\frac{k_1+k_2}{m_2}x_2+\frac{b_1}{m_2}x_3-\frac{b_1+b_2}{m_2}x_4-\frac{1}{m_2}u+(\frac{m_2k_2-b_1b_2-b_2^2}{m_2^2})w \\
	\end{aligned}
	\right.
\]

written in matrix:

\[
	\mathbf{\dot x}=
	\begin{bmatrix}
		0                 & 0                     & 1                 & 0                     \\
		0                 & 0                     & 0                 & 1                     \\
		-\dfrac{k_1}{m_1} & \dfrac{k_1}{m_1}      & -\dfrac{b_1}{m_1} & \dfrac{b_1}{m_1}      \\
		\dfrac{k_1}{m_2}  & -\dfrac{k_1+k_2}{m_2} & \dfrac{b_1}{m_2}  & -\dfrac{b_1+b_2}{m_2}
	\end{bmatrix}
	\mathbf{x} +
	\begin{bmatrix}
		0 \\ 0 \\ \dfrac{1}{m_1} \\-\dfrac{1}{m_2}
	\end{bmatrix}
	u+
	\begin{bmatrix}
		0 \\ \dfrac{b_2}{m_2}\\\dfrac{b_1b_2}{m_1m_2} \\ \dfrac{m_2k_2-b_1b_2-b_2^2}{m_2^2}
	\end{bmatrix}
	w
\]

\[
	y=x_1-x_2=
	\begin{bmatrix}
		1 & -1 & 0 & 0 \\
	\end{bmatrix}
	\mathbf{x}
\]

So the four matrices of the system could be defined as state matrix \(\bm{A}\), input matrix \(\bm{B}_1\), when disturbance input \(w=0\), input matrix \(\bm{B}_2\), when actuated force input \(u=0\), output matrix \(\bm{C}\) and direct transmission matrix \(\bm{D}\):

\[
	\bm{A}=
	\begin{bmatrix}
		0                 & 0                     & 1                 & 0                     \\
		0                 & 0                     & 0                 & 1                     \\
		-\dfrac{k_1}{m_1} & \dfrac{k_1}{m_1}      & -\dfrac{b_1}{m_1} & \dfrac{b_1}{m_1}      \\
		\dfrac{k_1}{m_2}  & -\dfrac{k_1+k_2}{m_2} & \dfrac{b_1}{m_2}  & -\dfrac{b_1+b_2}{m_2}
	\end{bmatrix},
\]
\[
	\bm{B}_1=
	\begin{bmatrix}
		0 \\ 0 \\ \dfrac{1}{m_1} \\-\dfrac{1}{m_2}
	\end{bmatrix}
	,
	\bm{B}_2=
	\begin{bmatrix}
		0 \\ \dfrac{b_2}{m_2}\\\dfrac{b_1b_2}{m_1m_2} \\\dfrac{m_2k_2-b_1b_2-b_2^2}{m_2^2}
	\end{bmatrix}
\]

\[
	\bm{C}=
	\begin{bmatrix}
		1 & -1 & 0 & 0
	\end{bmatrix}
	, \bm{D}=\bm{0}
\]
And here we have the state-space representation: 
\[
	\left \{
	\begin{aligned}
		\mathbf{\dot x} & =\bm{A}\mathbf{x}+\bm{B}_1u+\bm{B}_2w \\
		y               & =\bm{C}\mathbf{x}+\bm{D}u                 \\
	\end{aligned}
	\right.
\]

\section{System Solution in Matlab}

The three approaches to describe the dynamic mechanical system could be applied directly into Matlab, and the numerical and graphic solution of the system could be demonstrated with the help of the software.

\subsection{Transfer Function Equation}

According to the transfer function of the system,
\begin{align*}
	m_1 \ddot x_1 & =-b_1(\dot x_1 -\dot x_2 )-k_1(x_1-x_2)+u                                   \\
	m_2 \ddot x_2 & =b_1(\dot x_1 -\dot x_2 )+k_1(x_1-x_2)+b_2(\dot w - \dot x_2 )+k_2(w-x_2)-u
\end{align*}
and the actual value of the parameters, \(m_1=2500, m_2=320, k_1=80000, k_2=500000,b_1=350, b_2=15020\), both the input response under unit step actuated force and disturbance response under \SI{0.1}{\meter} disturbance could be determined.

\subsubsection{Input Response}
In order to obtain the input response of the system, we should input the following code into the command window of Matlab:

\begin{lstlisting}
	m1=2500; m2=320; k1=80000; k2=500000; b1 = 350; b2 = 15020;
	nump=[(m1+m2) b2 k2];
	denp=[(m1*m2) (m1*(b1+b2))+(m2*b1) (m1*(k1+k2))+(m2*k1)+(b1*b2) (b1*k2)+(b2*k1) k1*k2];
	'G(s)1'
	printsys(nump,denp)
	t=0:0.01:50;
	y= step(nump,denp,t);
	plot(t,y)
	grid
	title('open-loop response to unit step actuated force (transfer function)')
\end{lstlisting}


Then the software will respond by printing the transfer function of the system \(G_1(s)\) on the screen and showing a diagram of the systems response to the unit step actuated force as Figure\ \ref{fig:tfe_ir} shown.

\begin{lstlisting}
	ans =
	
		'G(s)1'
	
	 
	num/den = 
	 
							 2820 s^2 + 15020 s + 500000
	   -----------------------------------------------------------------------
	   800000 s^4 + 38537000 s^3 + 1480857000 s^2 + 1376600000 s + 40000000000
\end{lstlisting}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{matlabcode/tfe_ir.pdf}
	\caption{System Response to the Unit Step Actuated Force}\label{fig:tfe_ir}
\end{figure}

\subsubsection{Disturbance Response}

The system response to \(\SI{0.1}{\meter}\) step disturbance could also be calculated with the following code:

\begin{lstlisting}
	m1=2500; m2=320; k1=80000; k2=500000; b1 = 350; b2 = 15020;
	num2=[-(m1*b2) -(m1*k2) 0 0];
	den2=[(m1*m2) (m1*(b1+b2))+(m2*b1) (m1*(k1+k2))+(m2*k1)+(b1*b2) (b1*k2)+(b2*k1) k1*k2];
	'G(s)2'
	printsys(0.1*num2,den2)
	t=0:0.01:50;
	y2= step(0.1*num2,den2,t);
	plot(t,y2)
	grid
	title('open-loop response to 0.1 m step disturbance (transfer function)')
\end{lstlisting}

The software respond with the transfer function and the response curve as Figure\ \ref{fig:tfe_dr}.

\begin{lstlisting}
ans =

    'G(s)2'

 
num/den = 
 
                         -3755000 s^3 - 125000000 s^2
   -----------------------------------------------------------------------
   800000 s^4 + 38537000 s^3 + 1480857000 s^2 + 1376600000 s + 40000000000
\end{lstlisting}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{matlabcode/tfe_dr.pdf}
	\caption{System Response to Step Disturbance}\label{fig:tfe_dr}
\end{figure}



\subsection{State Space Matrix}
Based on the four matrices the system and the parameters, we could utilize the state space approach of Matlab to determine the two types of response of the system as well.


\subsubsection{Input Response}
When taking into account only the input \(u\) of the system, the input matrix \(\bm{B}_2\) of the disturbance is negligible. The code are as follows:

\begin{lstlisting}
	m1 = 2500; m2 = 320; k1 = 80000; k2 = 500000; b1 = 350; b2 = 15020;
	A = [
		0, 0, 1, 0;
		0, 0, 0, 1;
		-k1 / m1, k1 / m1, -b1 / m1, b1 / m1;
		k1 / m2, -(k1 + k2) / m2, b1 / m2, -(b1 + b2) / m2;
		];
	B = [0; 0; 1 / m1; -1 / m2];
	C = [1 -1 0 0];
	D = 0;
	sys = ss(A, B, C, D);
	t=0:0.01:50;
	y=step(sys,t);
	plot(t,y)
	grid
	title('open-loop response to unit step actuated force (state space)')	
\end{lstlisting}

The software responds directly with the Figure\ \ref{fig:ssm_ir}.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{matlabcode/ssm_ir.pdf}
	\caption{System Response to the Unit Step Actuated Force}\label{fig:ssm_ir}
\end{figure}

\subsubsection{Disturbance Response}
In order to determine the system response under \SI{0.1}{ \meter} step disturbance, we must assume that the input \(u=0\), thus the input matrix \(\bm{B}_1\) is ignored. The code is shown below:

\begin{lstlisting}
	m1 = 2500; m2 = 320; k1 = 80000; k2 = 500000; b1 = 350; b2 = 15020;
	A = [
		0, 0, 1, 0;
		0, 0, 0, 1;
		-k1 / m1, k1 / m1, -b1 / m1, b1 / m1;
		k1 / m2, -(k1 + k2) / m2, b1 / m2, -(b1 + b2) / m2;
		];
	B = [0; b2/m2; (b1*b2) /( m1*m2); (m2*k2-b1*b2-b2^2) / m2^2];
	C = [1 -1 0 0];
	D = 0;
	sys = ss(A, B, C, D);
	t=0:0.01:50;
	y=step(0.1*sys,t);
	plot(t,y)
	grid
	title('open-loop response to 0.1 m step disturbance (state space)')
\end{lstlisting}

The software executes the program and displays the disturbance response curve as Figure \ \ref{fig:ssm_dr} shows.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{matlabcode/ssm_dr.pdf}
	\caption{System Response to Step Disturbance}\label{fig:ssm_dr}
\end{figure}

\section{System Solution in Simulink}

With the help of Simulink in Matlab, the flow of the signal of the system could be demonstrated clearly in diagram connection.

\subsection{Deferential Equation}
Because of the complexity of the representation of the system in deferential equation, the connection of blocks in Simulink is much more complicated than that of the other two types of representation, namely the transfer function equation and state space matrix. All the parameters of the single blocks are set according to the parameters of the deferential equation and the known actual value of them. The deferential equation and block connection is shown in Figure\ \ref{fig:demodel}.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{demodel.pdf}
	\caption{the Block Connection Corresponding to the Deferential Equation}\label{fig:demodel}
\end{figure}

\subsubsection{Input Response}
The Step value of the disturbance generator is set to \(0\) in order to be ignored. The final output in the Scope block is shown in Figure\ \ref{fig:de_ir}.

\subsubsection{Disturbance Response}
At this time, all the parameters of input Step block is set to 0 for they don’t play any role in the disturbance response. The final output in the Scope block is shown in Figure\ \ref{fig:de_dr}.

\begin{figure}[htbp]
	\centering \subfigure[System Response to the Unit Step Actuated Force]{\includegraphics[width=0.48\textwidth]{matlabcode/de_ir.pdf}\label{fig:de_ir}}
	\centering \subfigure[System Response to Step Disturbance]{\includegraphics[width=0.48\textwidth]{matlabcode/de_dr.pdf}\label{fig:de_dr}}
	\caption{System Response with the Deferential Equation in Simulink}
\end{figure}

\subsection{Transfer Function Equation}
The system could be simplified into a signal generator, the transfer function of the system (input and disturbance transfer function respectively), and a scope, with whose help we could visualize the system response, and the connection of the blocks is shown in Figure\ \ref{fig:tfmodel}. The responses of input and disturbance are shown in Figure \ \ref{fig:tf_ir} and Figure \ \ref{fig:tf_dr}.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{tfmodel.pdf}
	\caption{the Block Connection Corresponding to the Transfer Function}\label{fig:tfmodel}
\end{figure}

\begin{figure}[htbp]
	\centering \subfigure[System Response to the Unit Step Actuated Force]{\includegraphics[width=0.48\textwidth]{matlabcode/de_ir.pdf}\label{fig:tf_ir}}
	\centering \subfigure[System Response to Step Disturbance]{\includegraphics[width=0.48\textwidth]{matlabcode/de_dr.pdf}\label{fig:tf_dr}}
	\caption{System Response with the Transfer Function in Simulink}
\end{figure}

\subsection{State Space Matrix}
The implementation of the state space approach in Simulink is the same as that of transfer function, expect the substitution of Block Transfer with block State-Space. The schematic connection of the system with state space approach is illustrated by the Figure\ \ref{fig:ssm}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{ssm.pdf}
	\caption{the Block Connection Corresponding to the State Space Matrix}\label{fig:ssm}
\end{figure}

\subsubsection{Input Response}

The settings of the block State-Space are shown in the following code block and we can see the output in Figure\ \ref{fig:simssm_ir}.
\begin{lstlisting}
	A = [
		0 0 1 0;
		0 0 0 1;
		-32 32 -0.14 0.14 ;
		250 -1812.5 1.09375 -48.03125
	]
	B = [0 ;0 ;0.0004 ;-0.003125]
	C = [1 -1 0 0]
	D = 0
\end{lstlisting}

\subsubsection{Disturbance Response}
The disturbance response is determined by the same procedure. The input matrix \(\bm{B}_1\) is replaced by the input matrix \(\bm{B}_2\) when we only pay attention to the disturbance response of the system. Here we demonstrate the configuration of \(\bm{B}_2\) and the final output are shown in Figure\ \ref{fig:simssm_dr}.

\begin{lstlisting}
	B = [0 ;46.9375; 6.57125; -691.966796875]
\end{lstlisting}

\begin{figure}[htbp]
	\centering \subfigure[System Response to the Unit Step Actuated Force]{\includegraphics[width=0.48\textwidth]{matlabcode/de_ir.pdf}\label{fig:simssm_ir}}
	\centering \subfigure[System Response to Step Disturbance]{\includegraphics[width=0.48\textwidth]{matlabcode/de_dr.pdf}\label{fig:simssm_dr}}
	\caption{System Response with the State Space Matrix in Simulink}
\end{figure}

\section{Design of a Feedback Controller}

From the graph of open-loop response for \SI{0.1}{\meter} step disturbance, we can see that when the bus passes a \SI{10}{\centi\meter} high bump on the road, the bus body will oscillate for an unacceptably long time (100 seconds) with larger amplitude, \SI{13}{\centi\meter}, than the initial impact. People sitting in the bus will not be comfortable with such an oscillation. The big overshoot (from the impact itself) and the slow settling time will cause damage to the suspension system. The solution to this problem is to add a feedback controller into the system to improve the performance.


\subsection{Schematic Diagram of the Closed Loop System}
According to the classic control theory, the performance of the system could be improved by adding an PID controller as closed loop feedback into the system. Since the amplitude of the disturbance response is about as 4000 times as the of input response, we added a Gain which is \(-10 000\) before the PID controller. Figure\ \ref{fig:pid} is the schematic diagram of the system.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{ssmmodel.pdf}
	\caption{the Block Connection of a Feedback Controller}\label{fig:pid}
\end{figure}

\subsection{Parameters of the PID Controller}
After several times of experiments, in which dozens of groups of all proportional, integral and derivative parameters are used in the PID Controller block, a final group of parameters of the PID Controller is determined. We configure the block as \(P=10, I=2, D=30\).


\subsection{Achieved System Performance under \SI{0.1}{\meter} Disturbance}
Using Simulink simulation, the final disturbance response of the system with the help of the added PID closed loop controller under \(\SI{0.1}{\meter}\) disturbance could be determined, the result is shown in Figure\ \ref{fig:result}.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.8\textwidth]{matlabcode/pid.pdf}
	\caption{Achieved System Performance}\label{fig:result}
\end{figure}

As we can see in the response curve under disturbance input, the system could achieve a better response than before. Here we denote the \emph{overshoot} as \( \sigma \) and \emph{settling time} as \(t_s\). The performance parameters of the ameliorated system are respectively:

\[ \sigma \approx 4.5 \% \]
\[ t_s \approx \SI{1.8}{s} \]

So both the overshoot and settling time could satisfy the requirements of the design task and the passengers will be comfortable even when the bus passes a \(\SI{10}{cm}\) high bump on the road. The suspension system of the vehicle could be effectively protected as well.

\end{document}
